//
// Created by Ivor_ on 2022/1/22.
//

#include "matrixOperate.h"

void prtCXY(_Complex double z0, double format) {
    char string[100] = "%";
    char nums[6];
    sprintf(nums, "%5.2f", format);
    strcat(string, nums);
    if (cimag(z0) > 0) {
        strcat(string, "f + %");
    } else {
        strcat(string, "f - %");
    }
    strcat(string, nums);
    strcat(string, "fi");
    printf(string, creal(z0), fabs(cimag(z0)));
}

void prtCRT(_Complex double z0, double format) {
    char string[100] = "%";
    char nums[6];
    sprintf(nums, "%5.2f", format);
    strcat(string, nums);
    strcat(string, "f * exp(%");
    strcat(string, nums);
    strcat(string, "fi)");
    printf(string, cabs(z0), carg(z0));
}

void fprtCXY(FILE* fp, _Complex double z0, double format) {
    char string[100] = "%";
    char nums[6];
    sprintf(nums, "%5.2f", format);
    strcat(string, nums);
    if (cimag(z0) > 0) {
        strcat(string, "f + %");
    } else {
        strcat(string, "f - %");
    }
    strcat(string, nums);
    strcat(string, "fi");
    fprintf(fp, string, creal(z0), fabs(cimag(z0)));
}

void fprtCRT(FILE* fp, _Complex double z0, double format) {
    char string[100] = "%";
    char nums[6];
    sprintf(nums, "%5.2f", format);
    strcat(string, nums);
    strcat(string, "f * exp(%");
    strcat(string, nums);
    strcat(string, "fi)");
    fprintf(fp, string, cabs(z0), carg(z0));
}

int* simDim(int len, ...) {
    int* dims = (int*)calloc(len, sizeof(int));
    va_list vl;
    va_start(vl, len);
    for (int i = 0; len; --len, ++i) {
        dims[i] = va_arg(vl, int);
    }
    va_end(vl);
    return dims;
}

double* simDVec(int dim, ...) {
    double* array = (double*)calloc(dim, sizeof(double));
    va_list vl;
    va_start(vl, dim);
    for (int i = 0; dim; --dim, ++i) {
        array[i] = va_arg(vl, double);
    }
    va_end(vl);
    return array;
}

_Complex double* simCVec(int dim, ...) {
    _Complex double* array = (_Complex double*)calloc(dim, sizeof(_Complex double));
    va_list vl;
    va_start(vl, dim);
    for (int i = 0; dim; --dim, ++i) {
        array[i] = va_arg(vl, _Complex double);
    }
    va_end(vl);
    return array;
}

double dLenVec(int dim, double* vec) {
    double len = 0.;
    for (int i = 0; i < dim; ++i) {
        len += pow(vec[i], 2.);
    }
    return len;
}

_Complex double cLenVec(int dim, _Complex double* vec) {
    double len = 0.;
    for (int i = 0; i < dim; ++i) {
        len += vec[i]* conj(vec[i]);
    }
    return sqrt(len) + 0.i;
}

double dInnVec(int dim, double* vec1, double* vec2) {
    double inn = 0.;
    for (int i = 0; i < dim; ++i) {
        inn += (vec1[i]*vec2[i]);
    }
    return inn;
}

_Complex double cInnVec(int dim, _Complex double* vec1, _Complex double* vec2) {
    _Complex double inn = 0. + 0.i;
    for (int i = 0; i < dim; ++i) {
        inn = inn + (conj(vec1[i])*vec2[i]);
    }
    return inn;
}

double** dOutVec(int* dims, double* vec1, double* vec2) {
    double** out = (double**)calloc(dims[0], sizeof(double*));
    for (int i = 0; i < dims[0]; ++i) {
        out[i] = (double*)calloc(dims[1], sizeof(double));
        for (int j = 0; j < dims[1]; ++j) {
            out[i][j] = (vec1[i]*vec2[j]);
        }
    }
    return out;
}

_Complex double** cOutVec(int* dims, _Complex double* vec1, _Complex double* vec2) {
    _Complex double** out = (_Complex double**)calloc(dims[0], sizeof(_Complex double*));
    for (int i = 0; i < dims[0]; ++i) {
        out[i] = (_Complex double*)calloc(dims[1], sizeof(_Complex double));
        for (int j = 0; j < dims[1]; ++j) {
            out[i][j] = vec1[i]*conj(vec2[j]);
        }
    }
    return out;
}

double* dStretchVec(int dim, double* vec, double factor) {
    double* steVec = (double*)calloc(dim, sizeof(_Complex double));
    for (int i = 0; i < dim; ++i) {
        steVec[i] = factor*vec[i];
    }
    return steVec;
}

_Complex double* cStretchVec(int dim, _Complex double* vec, _Complex double factor) {
    _Complex double* steVec = (_Complex double*)calloc(dim, sizeof(_Complex double));
    for (int i = 0; i < dim; ++i) {
        steVec[i] = factor*vec[i];
    }
    return steVec;
}

double* dGetZeroVec(int dim) {
    double* zeros = (double*)calloc(dim, sizeof(double));
    for (int i = 0; i < dim; ++i) {
        zeros[i] = 0.;
    }
    return zeros;
}

_Complex double* cGetZeroVec(int dim) {
    _Complex double* zeros = (_Complex double*)calloc(dim, sizeof(_Complex double));
    for (int i = 0; i < dim; ++i) {
        zeros[i] = 0. + 0.i;
    }
    return zeros;
}

double* dGetOnesVec(int dim) {
    double* ones = (double*)calloc(dim, sizeof(double));
    for (int i = 0; i < dim; ++i) {
        ones[i] = 1.;
    }
    return ones;
}

_Complex double* cGetOnesRVec(int dim) {
    _Complex double* onesR = (_Complex double*)calloc(dim, sizeof(_Complex double));
    for (int i = 0; i < dim; ++i) {
        onesR[i] = 1. + 0.i;
    }
    return onesR;
}

_Complex double* cGetOnesIVec(int dim) {
    _Complex double* onesI = (_Complex double*)calloc(dim, sizeof(_Complex double));
    for (int i = 0; i < dim; ++i) {
        onesI[i] = 0. + 1.i;
    }
    return onesI;
}

_Complex double* cGetOnesRIVec(int dim) {
    _Complex double* onesRI = (_Complex double*)calloc(dim, sizeof(_Complex double));
    for (int i = 0; i < dim; ++i) {
        onesRI[i] = 1. + 1.i;
    }
    return onesRI;
}

_Complex double* getCplxRVec(int dim, double* reVec) {
    return cNewV(2).new2(dim, reVec, dNewV(0).create(dim));
}

_Complex double* getCplxIVec(int dim, double* imVec) {
    return cNewV(2).new2(dim, dNewV(0).create(dim), imVec);
}

_Complex double* getCplxRIVec(int dim, double* reVec, double* imVec) {
    _Complex double* cVec = (_Complex double*)calloc(dim, sizeof(_Complex double));
    for (int i = 0; i < dim; ++i) {
        cVec[i] = reVec[i] + imVec[i]*1.i;
    }
    return cVec;
}

double* getVecRealPart(int dim, _Complex double* vec) {
    double* rv = (double*)calloc(dim, sizeof(double));
    for (int i = 0; i < dim; ++i) {
        rv[i] = creal(vec[i]);
    }
    return rv;
}

double* getVecImagPart(int dim, _Complex double* vec) {
    double* iv = (double*)calloc(dim, sizeof(double));
    for (int i = 0; i < dim; ++i) {
        iv[i] = cimag(vec[i]);
    }
    return iv;
}

double* getVecNormPart(int dim, _Complex double* vec) {
    double* mv = (double*)calloc(dim, sizeof(double));
    for (int i = 0; i < dim; ++i) {
        mv[i] = cabs(vec[i]);
    }
    return mv;
}

double* getVecArgPart(int dim, _Complex double* vec) {
    double* av = (double*)calloc(dim, sizeof(double));
    for (int i = 0; i < dim; ++i) {
        av[i] = carg(vec[i]);
    }
    return av;
}

double** dTrans(int* dims, double** mat) {
    double** trans = (double**)calloc(dims[1], sizeof(double*));
    for (int i = 0; i < dims[1]; ++i) {
        trans[i] = (double*)calloc(dims[0], sizeof(double));
        for (int j = 0; j < dims[0]; ++j) {
            trans[i][j] = mat[j][i];
        }
    }
    return trans;
}

_Complex double** cTrans(int* dims, _Complex double** mat) {
    _Complex double** trans = (_Complex double**)calloc(dims[1], sizeof(_Complex double*));
    for (int i = 0; i < dims[1]; ++i) {
        trans[i] = (_Complex double*)calloc(dims[0], sizeof(_Complex double));
        for (int j = 0; j < dims[0]; ++j) {
            trans[i][j] = mat[j][i];
        }
    }
    return trans;
}

struct dVecFeature dvf(void) {
    struct dVecFeature vf;
    vf.length = dLenVec;
    vf.innerProd = dInnVec;
    return vf;
}

struct cVecFeature cvf(void) {
    struct cVecFeature vf;
    vf.length = cLenVec;
    vf.innerProd = cInnVec;
    return vf;
}

struct dNewVec dNewV(int mode) {
    struct dNewVec dnv = {dGetZeroVec, getVecRealPart};
    if (((mode >= 0) && (mode <= 9)) || ((mode >= '0') && (mode <= '9'))) {
        switch (mode) {
            case 1:
            case '1':
                dnv.create = dGetOnesVec;
                break;
            default:
                dnv.create = dGetZeroVec;
                break;
        }
    } else {
        switch (mode) {
            case 'r':
            case 'R':
                dnv.sample = getVecRealPart;
                break;
            case 'i':
            case 'I':
                dnv.sample = getVecImagPart;
                break;
            case 'n':
            case 'N':
                dnv.sample = getVecNormPart;
                break;
            case 'a':
            case 'A':
                dnv.sample = getVecArgPart;
                break;
            default:
                fprintf(stderr, "No such mode in new vector, try to use R I N A or create method: %d\n", errno);
                exit(-1);
                break;
        }
    }
    return dnv;
}

struct cNewVec cNewV(int mode) {
    struct cNewVec cnv = {cGetZeroVec, getCplxRVec, getCplxRIVec};
    if (mode == 'n' || mode == 'N') {
        switch (mode) {
            case 'n':
                cnv.new1 = getCplxIVec;
                break;
            case 'N':
                cnv.new1 = getCplxRVec;
                break;
            default:
                fprintf(stderr, "No such mode in new vector, try to use n N or new2 method: %d\n", errno);
                exit(-1);
                break;
        }
    } else {
        switch (mode) {
            case 0:
            case '0':
                cnv.create = cGetZeroVec;
                break;
            case 1:
                cnv.create = cGetOnesRVec;
                break;
            case 'i':
                cnv.create = cGetOnesIVec;
                break;
            case '1':
                cnv.create = cGetOnesRIVec;
                break;
            default:
                cnv.new2 = getCplxRIVec;
                break;
        }
    }
    return cnv;
}

_Complex double** conjugateMat(int* dims, _Complex double** mat) {
    _Complex double** trans = (_Complex double**)calloc(dims[0], sizeof(_Complex double*));
    for (int i = 0; i < dims[0]; ++i) {
        trans[i] = (_Complex double*)calloc(dims[1], sizeof(_Complex double));
        for (int j = 0; j < dims[1]; ++j) {
            trans[i][j] = conj(mat[i][j]);
        }
    }
    return trans;
}

_Complex double** dagger(int* dims, _Complex double** mat) {
    return cCvtM().trans(dims, cCvtM().conj(dims, mat));
}

double** dGetZeroMat(int* dims) {
    double** zeros = (double**)calloc(dims[0], sizeof(double*));
    for (int i = 0; i < dims[0]; ++i) {
        zeros[i] = (double*)calloc(dims[1], sizeof(double));
        for (int j = 0; j < dims[1]; ++j) {
            zeros[i][j] = 0.;
        }
    }
    return zeros;
}

_Complex double** cGetZeroMat(int* dims) {
    _Complex double** zeros = (_Complex double**)calloc(dims[0], sizeof(_Complex double*));
    for (int i = 0; i < dims[0]; ++i) {
        zeros[i] = (_Complex double*)calloc(dims[1], sizeof(_Complex double));
        for (int j = 0; j < dims[1]; ++j) {
            zeros[i][j] = 0. + 0.i;
        }
    }
    return zeros;
}

double** dGetOnesMat(int* dims) {
    double** ones = (double**)calloc(dims[0], sizeof(double*));
    for (int i = 0; i < dims[0]; ++i) {
        ones[i] = (double*)calloc(dims[1], sizeof(double));
        for (int j = 0; j < dims[1]; ++j) {
            ones[i][j] = 1.;
        }
    }
    return ones;
}

_Complex double** cGetOnesRMat(int* dims) {
    _Complex double** onesR = (_Complex double**)calloc(dims[0], sizeof(_Complex double*));
    for (int i = 0; i < dims[0]; ++i) {
        onesR[i] = (_Complex double*)calloc(dims[1], sizeof(_Complex double));
        for (int j = 0; j < dims[1]; ++j) {
            onesR[i][j] = 1. + 0.i;
        }
    }
    return onesR;
}

_Complex double** cGetOnesIMat(int* dims) {
    _Complex double** onesI = (_Complex double**)calloc(dims[0], sizeof(_Complex double*));
    for (int i = 0; i < dims[0]; ++i) {
        onesI[i] = (_Complex double*)calloc(dims[1], sizeof(_Complex double));
        for (int j = 0; j < dims[1]; ++j) {
            onesI[i][j] = 0. + 1.i;
        }
    }
    return onesI;
}

_Complex double** cGetOnesRIMat(int* dims) {
    _Complex double** onesRI = (_Complex double**)calloc(dims[0], sizeof(_Complex double*));
    for (int i = 0; i < dims[0]; ++i) {
        onesRI[i] = (_Complex double*)calloc(dims[1], sizeof(_Complex double));
        for (int j = 0; j < dims[1]; ++j) {
            onesRI[i][j] = 1. + 1.i;
        }
    }
    return onesRI;
}

double** dGetUniteMat(int dim) {
    double** ones = (double**)calloc(dim, sizeof(double*));
    for (int i = 0; i < dim; ++i) {
        ones[i] = (double*)calloc(dim, sizeof(double));
        ones[i][i] = 1.;
    }
    return ones;
}

_Complex double** cGetUniteRMat(int dim) {
    _Complex double** unite = (_Complex double**)calloc(dim, sizeof(_Complex double*));
    for (int i = 0; i < dim; ++i) {
        unite[i] = (_Complex double*)calloc(dim, sizeof(_Complex double));
        unite[i][i] = 1. + 0.i;
    }
    return unite;
}

_Complex double** cGetUniteIMat(int dim) {
    _Complex double** unite = (_Complex double**)calloc(dim, sizeof(_Complex double*));
    for (int i = 0; i < dim; ++i) {
        unite[i] = (_Complex double*)calloc(dim, sizeof(_Complex double));
        unite[i][i] = 0. + 1.i;
    }
    return unite;
}

_Complex double** cGetUniteRIMat(int dim) {
    _Complex double** unite = (_Complex double**)calloc(dim, sizeof(_Complex double*));
    for (int i = 0; i < dim; ++i) {
        unite[i] = (_Complex double*)calloc(dim, sizeof(_Complex double));
        unite[i][i] = 1. + 1.i;
    }
    return unite;
}

_Complex double** getCplxRMat(int* dims, double** reMat) {
    return cNewM(SHAPE_RECTANGLE, 0).new0(dims, reMat, dNewM(SHAPE_RECTANGLE, 0).create_rm(dims));
}

_Complex double** getCplxIMat(int* dims, double** imMat) {
    return cNewM(SHAPE_RECTANGLE, 0).new0(dims, dNewM(SHAPE_RECTANGLE, 0).create_rm(dims), imMat);
}

_Complex double** getCplxRIMat(int* dims, double** reMat, double** imMat) {
    _Complex double** cMat = (_Complex double**)calloc(dims[0], sizeof(_Complex double*));
    for (int i = 0; i < dims[0]; ++i) {
        cMat[i] = (_Complex double*)calloc(dims[1], sizeof(_Complex double));
        for (int j = 0; j < dims[j]; ++j) {
            cMat[i][j] = reMat[i][j] + imMat[i][j]*1.i;
        }
    }
    return cMat;
}

double** getMatRealPart(int* dims, _Complex double** mat) {
    double** rm = (double**)calloc(dims[0], sizeof(double*));
    for (int i = 0; i < dims[0]; ++i) {
        rm[i] = (double*)calloc(dims[1], sizeof(double));
        for (int j = 0; j < dims[1]; ++j) {
            rm[i][j] = creal(mat[i][j]);
        }
    }
    return rm;
}

double** getMatImagPart(int* dims, _Complex double** mat) {
    double** im = (double**)calloc(dims[0], sizeof(double*));
    for (int i = 0; i < dims[0]; ++i) {
        im[i] = (double*)calloc(dims[1], sizeof(double));
        for (int j = 0; j < dims[1]; ++j) {
            im[i][j] = cimag(mat[i][j]);
        }
    }
    return im;
}

double** getMatNormPart(int* dims, _Complex double** mat) {
    double** mm = (double**)calloc(dims[0], sizeof(double*));
    for (int i = 0; i < dims[0]; ++i) {
        mm[i] = (double*)calloc(dims[1], sizeof(double));
        for (int j = 0; j < dims[1]; ++j) {
            mm[i][j] = cabs(mat[i][j]);
        }
    }
    return mm;
}

double** getMatArgPart(int* dims, _Complex double** mat) {
    double** am = (double**)calloc(dims[0], sizeof(double*));
    for (int i = 0; i < dims[0]; ++i) {
        am[i] = (double*)calloc(dims[1], sizeof(double));
        for (int j = 0; j < dims[1]; ++j) {
            am[i][j] = carg(mat[i][j]);
        }
    }
    return am;
}

double** dDiagMat(int dim, double* values) {
    double** mat = (double**)calloc(dim, sizeof(double*));
    for (int i = 0; i < dim; ++i) {
        mat[i] = (double*)calloc(dim, sizeof(double));
        mat[i][i] = values[i];
    }
    return mat;
}

_Complex double** cDiagMat(int dim, _Complex double* values) {
    _Complex double** mat = (_Complex double**)calloc(dim, sizeof(_Complex double*));
    for (int i = 0; i < dim; ++i) {
        mat[i] = (_Complex double*)calloc(dim, sizeof(_Complex double));
        mat[i][i] = values[i];
    }
    return mat;
}

double* * dVecCpy(int dim, double* vector) {
    double* * cpv = (double**)malloc(sizeof(double*));
    *cpv = (double*)calloc(dim, sizeof(double));
    for (int i = 0; i < dim; ++i) {
        (*cpv)[i] = vector[i];
    }
    return cpv;
}

_Complex double* * cVecCpy(int dim, _Complex double* vector) {
    _Complex double* * cpv = (_Complex double**)malloc(sizeof(_Complex double*));
    *cpv = (_Complex double*)calloc(dim, sizeof(_Complex double));
    for (int i = 0; i < dim; ++i) {
        (*cpv)[i] = vector[i];
    }
    return cpv;
}

double** * dMatCpy(int* dims, double** matrix) {
    double** * cpm = (double***)malloc(sizeof(double**));
    *cpm = (double**)calloc(dims[0], sizeof(double*));
    for (int i = 0; i < dims[0]; ++i) {
        (*cpm)[i] = (double*)calloc(dims[1], sizeof(double));
        for (int j = 0; j < dims[1]; ++j) {
            (*cpm)[i][j] = matrix[i][j];
        }
    }
    return cpm;
}

_Complex double** * cMatCpy(int* dims, _Complex double** matrix) {
    _Complex double** * cpm = (_Complex double***)malloc(sizeof(_Complex double**));
    *cpm = (_Complex double**)calloc(dims[0], sizeof(_Complex double*));
    for (int i = 0; i < dims[0]; ++i) {
        (*cpm)[i] = (_Complex double*)calloc(dims[1], sizeof(_Complex double));
        for (int j = 0; j < dims[1]; ++j) {
            (*cpm)[i][j] = matrix[i][j];
        }
    }
    return cpm;
}

struct dNewMat dNewM(int shape, int mode) {
    struct dNewMat dnm = {NULL, NULL, NULL, dDiagMat};
    if (shape == SHAPE_SQUARE) {
        switch (mode) {
            case 'U':
            case 'u':
                dnm.create_sm = dGetUniteMat;
                break;
            default:
                fprintf(stderr, "No such mode in new matrix_s, try to use U or matrix_r method: %d\n", errno);
                exit(-1);
                break;
        }
    } else if (shape == SHAPE_RECTANGLE) {
        switch (mode) {
            case 1:
            case '1':
                dnm.create_rm = dGetOnesMat;
                break;
            case 0:
            case '0':
                dnm.create_rm = dGetZeroMat;
                break;
            case 'R':
            case 'r':
                dnm.sample = getMatRealPart;
                break;
            case 'I':
            case 'i':
                dnm.sample = getMatImagPart;
                break;
            case 'N':
            case 'n':
                dnm.sample = getMatNormPart;
                break;
            case 'A':
            case 'a':
                dnm.sample = getMatArgPart;
                break;
            default:
                fprintf(stderr, "No such mode in new matrix_r, try to use 0 1 R I N A or matrix_s method: %d\n", errno);
                exit(-1);
                break;
        }
    } else {
        fprintf(stderr, "No such shape in matrix, try to input SHAPE_ and choose one: %d\n", errno);
        exit(-1);
    }
    return dnm;
}

struct cNewMat cNewM(int shape, int mode) {
    struct cNewMat cnm = {NULL, NULL, getCplxRIMat, NULL, NULL};
    if (shape == SHAPE_SQUARE) {
        switch (mode) {
            case 'U':
                cnm.create_sm = cGetUniteRMat;
                break;
            case 'u':
                cnm.create_sm = cGetUniteIMat;
                break;
            case 'd':
                cnm.create_sm = cGetUniteRIMat;
                break;
            case 'D':
                cnm.new2 = cDiagMat;
                break;
            default:
                fprintf(stderr, "No such mode in new matrix, try to use U u d D or matrix_r method: %d\n", errno);
                exit(-1);
                break;
        }
    } else if (shape == SHAPE_RECTANGLE) {
        switch (mode) {
            case 0:
            case '0':
                cnm.create_rm = cGetZeroMat;
                break;
            case 1:
                cnm.create_rm = cGetOnesRMat;
                break;
            case 'I':
            case 'i':
                cnm.create_rm = cGetOnesIMat;
                break;
            case '1':
                cnm.create_rm = cGetOnesRIMat;
                break;
            case 'R':
                cnm.new1 = getCplxRMat;
                break;
            case 'r':
                cnm.new1 = getCplxIMat;
                break;
            default:
                fprintf(stderr, "No such mode in new matrix, try to use 0 1 i R r or matrix_s method: %d\n", errno);
                exit(-1);
                break;
        }
    } else {
        fprintf(stderr, "No such shape in matrix, try to input SHAPE_ and choose one: %d\n", errno);
        exit(-1);
    }
    return cnm;
}

struct dConvertVec dCvtV(void) {
    struct dConvertVec dcv;
    dcv.stretch = dStretchVec;
    return dcv;
}

struct cConvertVec cCvtV(void) {
    struct cConvertVec ccv;
    ccv.stretch = cStretchVec;
    return ccv;
}

void dPrtVec(int dim, double* vec, double format) {
    for (int i = 0; i < dim; ++i) {
        char string[100] = "%";
        char nums[6];
        sprintf(nums, "%5.2f", format);
        strcat(string, nums);
        strcat(string, "f\t");
        printf(string, vec[i]);
    }
    putchar(10);
}

void dFprtVec(FILE* fp, int dim, double* vec, double format) {
    for (int i = 0; i < dim; ++i) {
        char string[100] = "%";
        char nums[6];
        sprintf(nums, "%5.2f", format);
        strcat(string, nums);
        strcat(string, "f\t");
        fprintf(fp, string, vec[i]);
    }
    fputc(10, fp);
}

void dPrtMat(int* dims, double** mat, double format) {
    for (int i = 0; i < dims[0]; ++i) {
        for (int j = 0; j < dims[1]; ++j) {
            char string[100] = "%";
            char nums[6];
            sprintf(nums, "%5.2f", format);
            strcat(string, nums);
            strcat(string, "f\t");
            printf(string, mat[i][j]);
        }
        putchar(10);
    }
}

void dFprtMat(FILE* fp, int* dims, double** mat, double format) {
    for (int i = 0; i < dims[0]; ++i) {
        for (int j = 0; j < dims[1]; ++j) {
            char string[100] = "%";
            char nums[6];
            sprintf(nums, "%5.2f", format);
            strcat(string, nums);
            strcat(string, "f\t");
            fprintf(fp, string, mat[i][j]);
        }
        fputc(10, fp);
    }
}

void cPrtVecXY(int dim, _Complex double* vec, double format) {
    for (int i = 0; i < dim; ++i) {
        prtCXY(vec[i], format);
        putchar('\t');
    }
    putchar(10);
}

void cFprtVecXY(FILE* fp, int dim, _Complex double* vec, double format) {
    for (int i = 0; i < dim; ++i) {
        fprtCXY(fp, vec[i], format);
        fputc('\t', fp);
    }
    fputc(10, fp);
}

void cPrtMatXY(int* dims, _Complex double** mat, double format) {
    for (int i = 0; i < dims[0]; ++i) {
        for (int j = 0; j < dims[1]; ++j) {
            prtCXY(mat[i][j], format);
            putchar('\t');
        }
        putchar(10);
    }
}

void cFprtMatXY(FILE* fp, int* dims, _Complex double** mat, double format) {
    for (int i = 0; i < dims[0]; ++i) {
        for (int j = 0; j < dims[1]; ++j) {
            fprtCXY(fp, mat[i][j], format);
            fputc('\t', fp);
        }
        fputc(10, fp);
    }
}

void cPrtVecRT(int dim, _Complex double* vec, double format) {
    for (int i = 0; i < dim; ++i) {
        prtCRT(vec[i], format);
        putchar('\t');
    }
    putchar(10);
}

void cFprtVecRT(FILE* fp, int dim, _Complex double* vec, double format) {
    for (int i = 0; i < dim; ++i) {
        fprtCRT(fp, vec[i], format);
        fputc('\t', fp);
    }
    fputc(10, fp);
}

void cPrtMatRT(int* dims, _Complex double** mat, double format) {
    for (int i = 0; i < dims[0]; ++i) {
        for (int j = 0; j < dims[1]; ++j) {
            prtCRT(mat[i][j], format);
            putchar('\t');
        }
        putchar(10);
    }
}

void cFprtMatRT(FILE* fp, int* dims, _Complex double** mat, double format) {
    for (int i = 0; i < dims[0]; ++i) {
        for (int j = 0; j < dims[1]; ++j) {
            fprtCRT(fp, mat[i][j], format);
            fputc('\t', fp);
        }
        fputc(10, fp);
    }
}

struct printArray prtArr(int form) {
    struct printArray pa;
    pa.dv = dPrtVec;
    pa.dm = dPrtMat;
    pa.fdv = dFprtVec;
    pa.fdm = dFprtMat;
    if (form == FORM_XY) {
        pa.cv = cPrtVecXY;
        pa.cm = cPrtMatXY;
        pa.fcv = cFprtVecXY;
        pa.fcm = cFprtMatXY;
    } else if (form == FORM_RT) {
        pa.cv = cPrtVecRT;
        pa.cm = cPrtMatRT;
        pa.fcv = cFprtVecRT;
        pa.fcm = cFprtMatRT;
    } else {
        fprintf(stderr, "No such form in print array, try to choose from FORM_ or input 0 for double type: %d\n", errno);
        exit(-1);
    }
    return pa;
}

void dExchangeLines(double** mat, int* dims, int* locates) {
    for (int i = 0; i < dims[1]; ++i) {
        double temp = mat[locates[0]][i];
        mat[locates[0]][i] = mat[locates[1]][i];
        mat[locates[1]][i] = temp;
    }
}

void dExchangeRows(double** mat, int* dims, int* locates) {
    for (int i = 0; i < dims[0]; ++i) {
        double temp = mat[i][locates[0]];
        mat[i][locates[0]] = mat[i][locates[1]];
        mat[i][locates[1]] = temp;
    }
}

void dMultipleLine(double** mat, int* dims, int locate, double factor) {
    for (int i = 0; i < dims[1]; ++i) {
        mat[locate][i] *= factor;
    }
}

void dMultipleRow(double** mat, int* dims, int locate, double factor) {
    for (int i = 0; i < dims[0]; ++i) {
        mat[i][locate] *= factor;
    }
}

void dDoubledLines(double** mat, int* dims, int* locates, double factor) {
    for (int i = 0; i < dims[1]; ++i) {
        mat[locates[1]][i] += (factor*mat[locates[0]][i]);
    }
}

void dDoubledRows(double** mat, int* dims, int* locates, double factor) {
    for (int i = 0; i < dims[0]; ++i) {
        mat[i][locates[1]] += (factor*mat[i][locates[0]]);
    }
}

void cExchangeLines(_Complex double** mat, int* dims, int* locates) {
    for (int i = 0; i < dims[1]; ++i) {
        _Complex double temp = mat[locates[0]][i];
        mat[locates[0]][i] = mat[locates[1]][i];
        mat[locates[1]][i] = temp;
    }
}

void cExchangeRows(_Complex double** mat, int* dims, int* locates) {
    for (int i = 0; i < dims[0]; ++i) {
        _Complex double temp = mat[i][locates[0]];
        mat[i][locates[0]] = mat[i][locates[1]];
        mat[i][locates[1]] = temp;
    }
}

void cMultipleLine(_Complex double** mat, int* dims, int locate, _Complex double factor) {
    for (int i = 0; i < dims[1]; ++i) {
        mat[locate][i] = mat[locate][i]*factor;
    }
}

void cMultipleRow(_Complex double** mat, int* dims, int locate, _Complex double factor) {
    for (int i = 0; i < dims[0]; ++i) {
        mat[i][locate] = mat[i][locate]*factor;
    }
}

void cDoubledLines(_Complex double** mat, int* dims, int* locates, _Complex double factor) {
    for (int i = 0; i < dims[1]; ++i) {
        mat[locates[1]][i] = mat[locates[1]][i] + (factor*mat[locates[0]][i]);
    }
}

void cDoubledRows(_Complex double** mat, int* dims, int* locates, _Complex double factor) {
    for (int i = 0; i < dims[0]; ++i) {
        mat[i][locates[1]] = mat[i][locates[1]] + (factor*mat[i][locates[0]]);
    }
}

struct dSwitchMat dSwtM(char mode) {
    struct dSwitchMat dsm;
    if (mode == 'l' || mode == 'L') {
        dsm.exc = dExchangeLines;
        dsm.mul = dMultipleLine;
        dsm.dbl = dDoubledLines;
    } else if (mode == 'r' || mode == 'R') {
        dsm.exc = dExchangeRows;
        dsm.mul = dMultipleRow;
        dsm.dbl = dDoubledRows;
    } else {
        fprintf(stderr, "No such form in convert matrix, try to use r or l mode: %d\n", errno);
        exit(-1);
    }
    return dsm;
}

struct cSwitchMat cSwtM(char mode) {
    struct cSwitchMat csm;
    if (mode == 'l' || mode == 'L') {
        csm.exc = cExchangeLines;
        csm.mul = cMultipleLine;
        csm.dbl = cDoubledLines;
    } else if (mode == 'r' || mode == 'R') {
        csm.exc = cExchangeRows;
        csm.mul = cMultipleRow;
        csm.dbl = cDoubledRows;
    } else {
        fprintf(stderr, "No such form in convert matrix, try to use r or l mode: %d\n", errno);
        exit(-1);
    }
    return csm;
}

struct dConvertMat dCvtM(void) {
    struct dConvertMat dcm;
    dcm.trans = dTrans;
    return dcm;
}

struct cConvertMat cCvtM(void) {
    struct cConvertMat ccm;
    ccm.trans = cTrans;
    ccm.conj = conjugateMat;
    ccm.dag = dagger;
    return ccm;
}

double** dMatrixLineSimAlpha(int* dims, double** matrix) {
    int minDim = dims[0] < dims[1] ? dims[0] : dims[1];
    double** * res = dMatCpy(dims, matrix);
    int main = -1;
    for (int j = 0; j < minDim; j++) {
        int i, flag = 0;
        for (i = j; i < minDim; i++) {
            if (fabs((*res)[i][j]) > 1.e-9) {
                flag = 1;
                main++;
                break;
            }
        }
        if (flag == 1) {
            if (i != main) {
                dSwtM('L').exc(*res, dims, simDim(2, i, main));
            }
            for (int k = 0; k < dims[0]; k++) {
                if (k != main) {
                    dSwtM('L').dbl(*res, dims, simDim(2, main, k), - (*res)[k][j]/(*res)[main][j]);
                }
            }
        }
    }
    return *res;
}

double** dMatrixLineSimBeta(int* dims, double** matrix) {
    double** res = dMatrixLineSimAlpha(dims, matrix);
    for (int i = 0; i < dims[0]; i++) {
        double number = 0.;
        for (int j = 0; j < dims[1]; j++) {
            if (fabs(number) > 1.e-9) {
                res[i][j] /= number;
            }
            else if (fabs(res[i][j]) > 1.e-9) {
                number = res[i][j];
                res[i][j] = 1.;
            }
        }
    }
    return res;
}

_Complex double** cMatrixLineSimAlpha(int* dims, _Complex double** matrix) {
    int minDim = dims[0] < dims[1] ? dims[0] : dims[1];
    _Complex double** * res = cMatCpy(dims, matrix);
    int main = -1;
    for (int j = 0; j < minDim; j++) {
        int i, flag = 0;
        for (i = j; i < minDim; i++) {
            if (cabs((*res)[i][j]) > 1.e-9) {
                flag = 1;
                main++;
                break;
            }
        }
        if (flag == 1) {
            if (i != main) {
                cSwtM('L').exc(*res, dims, simDim(2, i, main));
            }
            for (int k = 0; k < dims[0]; k++) {
                if (k != main) {
                    cSwtM('L').dbl(*res, dims, simDim(2, main, k), - (*res)[k][j]/(*res)[main][j]);
                }
            }
        }
    }
    return *res;
}

_Complex double** cMatrixLineSimBeta(int* dims, _Complex double** matrix) {
    _Complex double** res = cMatrixLineSimAlpha(dims, matrix);
    for (int i = 0; i < dims[0]; i++) {
        _Complex double number = 0.;
        for (int j = 0; j < dims[1]; j++) {
            if (cabs(number) > 1.e-9) {
                res[i][j] /= number;
            }
            else if (cabs(res[i][j]) > 1.e-9) {
                number = res[i][j];
                res[i][j] = 1.;
            }
        }
    }
    return res;
}

double** dMatrixRowSimAlpha(int* dims, double** matrix) {
    return dCvtM().trans(simDim(2, dims[1], dims[0]), dMatrixLineSimAlpha(simDim(2, dims[1], dims[0]), dCvtM().trans(dims, matrix)));
}

double** dMatrixRowSimBeta(int* dims, double** matrix) {
    return dCvtM().trans(simDim(2, dims[1], dims[0]), dMatrixLineSimBeta(simDim(2, dims[1], dims[0]), dCvtM().trans(dims, matrix)));
}

_Complex double** cMatrixRowSimAlpha(int* dims, _Complex double** matrix) {
    return cCvtM().trans(simDim(2, dims[1], dims[0]), cMatrixLineSimAlpha(simDim(2, dims[1], dims[0]), cCvtM().trans(dims, matrix)));
}

_Complex double** cMatrixRowSimBeta(int* dims, _Complex double** matrix) {
    return cCvtM().trans(simDim(2, dims[1], dims[0]), cMatrixLineSimBeta(simDim(2, dims[1], dims[0]), cCvtM().trans(dims, matrix)));
}

double** * dCreateGivensMatrix(int dim, int* locateIndex, double rotateAngle) {
    double** * retMat = (double***)malloc(sizeof(double**));
    *retMat = dNewM(SHAPE_SQUARE, 'U').create_sm(dim);
    (*retMat)[locateIndex[0]][locateIndex[0]] = cos(rotateAngle);
    (*retMat)[locateIndex[1]][locateIndex[1]] = cos(rotateAngle);
    (*retMat)[locateIndex[1]][locateIndex[0]] = -sin(rotateAngle);
    (*retMat)[locateIndex[0]][locateIndex[1]] = sin(rotateAngle);
    return retMat;
}

_Complex double** * cCreateGivensMatrix(int dim, int* locateIndex, double rotateAngle, double* phases) {
    _Complex double** * retMat = (_Complex double***)malloc(sizeof(_Complex double**));
    *retMat = cNewM(SHAPE_SQUARE, 'U').create_sm(dim);
    (*retMat)[locateIndex[0]][locateIndex[0]] = cos(rotateAngle)*cexp(1.i*phases[0]);
    (*retMat)[locateIndex[1]][locateIndex[1]] = cos(rotateAngle)*cexp(-1.i*phases[0]);
    (*retMat)[locateIndex[1]][locateIndex[0]] = -sin(rotateAngle)*cexp(-1.i*phases[1]);
    (*retMat)[locateIndex[0]][locateIndex[1]] = sin(rotateAngle)*cexp(1.i*phases[1]);
    return retMat;
}

void dGivensTransDirect(int dim, double** * matrixPointer, double** givensMat) {
    double** * tempMat = dMulN().multiple(simDim(4, dim, dim, dim, dim), *dMulN().inner(simDim(4, dim, dim, dim, dim), givensMat, *matrixPointer), givensMat);
    freeDMat(dim, *matrixPointer);
    *matrixPointer = *tempMat;
}

void cGivensTransDirect(int dim, _Complex double** * matrixPointer, _Complex double** givensMat) {
    _Complex double** * tempMat = cMulN().multiple(simDim(4, dim, dim, dim, dim), *cMulN().inner(simDim(4, dim, dim, dim, dim), givensMat, *matrixPointer), givensMat);
    freeCMat(dim, *matrixPointer);
    *matrixPointer = *tempMat;
}

void dGivensTrans(int dim, double** * matrixPointer, int* locateIndex, double rotateAngle) {
    double** * givensMat = (double***)malloc(sizeof(double**));
    *givensMat = dNewM(SHAPE_SQUARE, 'U').create_sm(dim);
    (*givensMat)[locateIndex[0]][locateIndex[0]] = cos(rotateAngle);
    (*givensMat)[locateIndex[1]][locateIndex[1]] = cos(rotateAngle);
    (*givensMat)[locateIndex[1]][locateIndex[0]] = -sin(rotateAngle);
    (*givensMat)[locateIndex[0]][locateIndex[1]] = sin(rotateAngle);
    double** * tempMat = dMulN().multiple(simDim(4, dim, dim, dim, dim), *dMulN().inner(simDim(4, dim, dim, dim, dim), *givensMat, *matrixPointer), *givensMat);
    freeDMat(dim, *givensMat);
    freeDMat(dim, *matrixPointer);
    *matrixPointer = *tempMat;
}

void cGivensTrans(int dim, _Complex double** * matrixPointer, int* locateIndex, double rotateAngle, double* phases) {
    _Complex double** * givensMat = (_Complex double***)malloc(sizeof(_Complex double**));
    *givensMat = cNewM(SHAPE_SQUARE, 'U').create_sm(dim);
    (*givensMat)[locateIndex[0]][locateIndex[0]] = cos(rotateAngle) * cexp(1.i * phases[0]);
    (*givensMat)[locateIndex[1]][locateIndex[1]] = cos(rotateAngle) * cexp(-1.i * phases[0]);
    (*givensMat)[locateIndex[1]][locateIndex[0]] = -sin(rotateAngle) * cexp(-1.i * phases[1]);
    (*givensMat)[locateIndex[0]][locateIndex[1]] = sin(rotateAngle) * cexp(1.i * phases[1]);
    _Complex double** * tempMat = cMulN().multiple(simDim(4, dim, dim, dim, dim), *cMulN().inner(simDim(4, dim, dim, dim, dim), *givensMat, *matrixPointer), *givensMat);
    freeCMat(dim, *givensMat);
    freeCMat(dim, *matrixPointer);
    *matrixPointer = *tempMat;
}

double** * dCreateHouseholderMatrix(int dim, double* vector) {
    double** * retMat = (double***)malloc(sizeof(double**));
    *retMat = dNewM(SHAPE_SQUARE, 'U').create_sm(dim);
    for (int i = 0; i < dim; ++i) {
        for (int j = 0; j < dim; ++j) {
            (*retMat)[i][j] -= 2.*vector[i]*vector[j];
        }
    }
    return retMat;
}

_Complex double** * cCreateHouseholderMatrix(int dim, _Complex double* vector) {
    _Complex double** * retMat = (_Complex double***)malloc(sizeof(_Complex double**));
    *retMat = cNewM(SHAPE_SQUARE, 'U').create_sm(dim);
    for (int i = 0; i < dim; ++i) {
        for (int j = 0; j < dim; ++j) {
            (*retMat)[i][j] -= 2.*vector[i]*conj(vector[j]);
        }
    }
    return retMat;
}

double** * dHouseholderReflex(int dim, double* vector, double* direction) {
    double* * normalizedVector = dVecCpy(dim, vector);
    double* * normalizedDirection = dVecCpy(dim, direction);
    double vecLen = dLenVec(dim, *normalizedVector);
    double dirLen = dLenVec(dim, *normalizedDirection);
    for (int i = 0; i < dim; ++i) {
        (*normalizedVector)[i] /= vecLen;
        (*normalizedDirection)[i] /= dirLen;
    }
    double* chv = dNewV(0).create(dim);
    for (int i = 0; i < dim; ++i) {
        chv[i] = (*normalizedVector)[i] - (*normalizedDirection)[i];
    }
    freeDVec(*normalizedVector);
    freeDVec(*normalizedDirection);
    double chvLen = dLenVec(dim, chv);
    for (int i = 0; i < dim; ++i) {
        chv[i] /= chvLen;
    }
    double** * retMat = dCreateHouseholderMatrix(dim, chv);
    freeDVec(chv);
    return retMat;
}

_Complex double** * cHouseholderReflex(int dim, _Complex double* vector, _Complex double* direction) {
    _Complex double* * normalizedVector = cVecCpy(dim, vector);
    _Complex double* * normalizedDirection = cVecCpy(dim, direction);
    _Complex double vecLen = cLenVec(dim, *normalizedVector);
    _Complex double dirLen = cLenVec(dim, *normalizedDirection);
    for (int i = 0; i < dim; ++i) {
        (*normalizedVector)[i] /= vecLen;
        (*normalizedDirection)[i] /= dirLen;
    }
    _Complex double inp = cInnVec(dim, *normalizedVector, *normalizedDirection);
    _Complex double* chv = cNewV(0).create(dim);
    for (int i = 0; i < dim; ++i) {
        chv[i] = (*normalizedVector)[i]*cexp(carg(inp)*1.i) - (*normalizedDirection)[i];
    }
    freeCVec(*normalizedVector);
    freeCVec(*normalizedDirection);
    _Complex double chvLen = cLenVec(dim, chv);
    _Complex double** * retMat;
    if (cabs(chvLen)/dim < 1.e-9) {
        retMat = (_Complex double***)malloc(sizeof(double**));
        *retMat = cNewM(SHAPE_SQUARE, 'U').create_sm(dim);
    } else {
        for (int i = 0; i < dim; ++i) {
            chv[i] /= chvLen;
        }
        retMat = cCreateHouseholderMatrix(dim, chv);
    }
    freeCVec(chv);
    for (int i = 0; i < dim; ++i) {
        for (int j = 0; j < dim; ++j) {
            (*retMat)[i][j] *= cexp(carg(inp)*1.i);
        }
    }
    return retMat;
}
